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Abstract 

Numerical approximations of shock waves sometimes suffer from instabil- 
ities called carbuncles. Techniques for suppressing carbuncles are trial- 
and-error and lack in reliability and generality, partly because theoreti- 
cal knowledge about carbuncles is equally unsatisfactory. It is not known 
which numerical schemes are affected in which circumstances, what causes 
carbuncles to appear and whether carbuncles are purely numerical arti- 
facts or rather features of a continuum equation or model. This work 
presents evidence towards the latter: it is conjectured that carbuncles are 
a special class of non-physical entropy solutions. Using a new technique 
for triggering a single carbuncle, its structure is computed in detail in 
similarity coordinates. 



1 Introduction 



Numerical calculations of shock waves in inviscid compressible flow sometimes 
suffer from instabilities, which were christened carbuncles in |PI88j . Figure ^ 
presents an example: a shock wave in front of a forward-facing step, normally 
expected to be smooth, develops several small structures called "carbuncles" 
in the literature. Figure |2 shows these carbuncles in more detail. For other 
examples of carbuncle phenomena see |PI88| , |DMG03| or |Qui94| Figure 5] . 

There is no clear rule for predicting which numerical methods are likely to be 
affected or in which circumstances carbuncles might appear (a summary of var- 
ious proposals and explanations can be found in |DMG03| ). Figure [T] displays 
two characteristic situations: the upper group of carbuncles is located near a 
boundary where numerical boundary layers act; the lower group appears in a 
region where the shock is almost but not exactly parallel to the grid edges. Some 
authors single out the Roe scheme |Roe81| as particularly prone to carbuncles 
whereas others report bad experiences with the Godunov method |God59j . It 
appears that the theoretical consensus on carbuncles is as volatile as the car- 
buncles themselves. However, recently researcher have focused on instability of 
discrete plane shock waves as an ingredient for the genesis of carbuncles (see 

e.g. ibivic ioai gftxnnni iamm) ). 
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As a consequence of the theoretical confusion, there are no safe recipes for elim- 
inating carbuncles. Although numerical viscosity in various forms can suppress 
the carbuncle phenomenon, the exact amount required is not known. Since 
numerical viscosity smears discontinuities and decreases overall accuracy, it is 
desirable to have precise theoretical criteria for carbuncle avoidance. 

Due to the nature of carbuncles as phenomena without precise mathematical 
definition, it is not clear whether there is unanimous agreement in the scientific 
community on which phenomena can be classified as carbuncles. Numerical 
calculations of discontinuities in gas flow can suffer from many other problems, 
such as linear instability due to excessive time steps or poor choices for opera- 
tor discretization, nonlinear phenomena in the presence of large total variation, 
Kelvin-Helmholtz instabilities, expansion shocks etc. We hope the reader will 
agree that the quasi-continuum carbuncle patterns computed in standard (Fig- 
ure |3l and similarity (Figure Q) coordinates in this paper correspond to to the 
small discrete carbuncles shown in Figures ^ and [21 or in the references above. 

All results displayed in this paper were done for the 2D isentropic Eulcr equa- 
tions with 7 := 7/5; however, analogous results can be generated for other 7 
and for the nonisentropic Euler equations as well. In all diagrams, flow into the 
domains is Mach 3.0 horizontally from left to right. Unless otherwise noted, all 
pictures show the horizontal component of velocity. 



2 Carbuncles in standard coordinates 

In this section we describe how carbuncles can be triggered in a large variety of 
numerical schemes in standard coordinates (here, "standard" means coordinates 
(t,x) as compared to similarity coordinates (t, j).) We start with a steady 
vertical plane shock wave computation in a uniform Cartesian grid. The upper 
and lower domain boundaries are solid walls with slip boundary condition; fluxes 
across the left resp. right domain boundary are computed with pre- resp. post- 
shock state in "pseudo-cells" on the outside of the domain boundaries. In this 
setting, the exact plane shock wave is reproduced faithfully by discrete shocks. 

Now the initial and boundary data is modified as follows: we choose an edge in 
the center of the left boundary; in its pseudo-cell and in all cells in a horizontal 
one-cell-high "filament" from this edge to the shock, the horizontal velocity is 
set to zero. This filament immediately triggers a structure (see Figure [3J that 
grows steadily (and eventually reaches the domain boundaries). This struc- 
ture, somewhat reminiscent of a stingray, is very similar to the small discrete 
carbuncles shown in Figure |21 

Figurc|3has been computed with the Godunov scheme |God59j : similar carbun- 
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Figure 1. Curved shock with car- Figure 2. Closeup of carbuncles in Fig- 

buncles in front of a forward- facing ure^ 

step 




Solid wall Solid wall 

Figure 3. A stingray-like structure resembling carbuncles 
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cles can be observed in the local Lax-Friedrichs |SQ89| or the Osher-Solomon 
|OS82| schemes. Higher-order schemes are likewise affected (this is not particu- 
larly surprising because the creation of the carbuncle takes place in nonsmooth 
or even discontinuous regions where higher-order schemes revert to first order). 
Moreover, adding generous amounts of numerical viscosity anywhere except near 
the filament and the carbuncle tip does not destroy the carbuncle. This suggests 
that the tip is the unphysical part of the carbuncle (unless of course carbuncles 
are physically meaningful). 

The reader may suspect that a pronounced disturbance like the carbuncle in 
Figure is inevitable given the reduced mass and momentum transport in the 
filament. However, the filament is too thin to support this reasoning. When the 
grid is refined and the filament height (and hence the L 1 norm perturbation to 
the initial and boundary data) decreases, the carbuncle does not diminuish; it 
always grows until it pushes the confines of the domain. 



3 Carbuncles in similarity coordinates 

The new method for triggering carbuncles is most effective in standard coor- 
dinates where it can be successfully applied to many numerical methods. It is 
observed that the carbuncle pattern grows at a constant velocity and that its 
inner structure settles in after some time. This suggests that carbuncles are 
self-similar flow patterns. Moreover, the steadily growing carbuncles eventu- 
ally reach the boundary of any spatial domain. Therefore, for computing their 
structure in detail similarity coordinates (t, f ) are more appropriate. 

A major drawback is that similarity coordinates add dissipation because grid 
cells move along the shock and smoothen the solution by averaging. For most 
numerical schemes, this decreases carbuncles to a point where it cannot be 
said for sure whether carbuncle size is independent of filament size. More im- 
portantly, it is not possible to create carbuncles of sufficient size to compute 
significant detail. However, the Godunov scheme does produce robust and large 
carbuncles in our setting. 

Here we compute on a Cartesian grid with adaptive refinement to achieve better 
resolution. Similarity coordinates are treated by the moving-edge modifications 
discussed in [EllOQI |E1105| . To save time, only the upper halfplane of the obvi- 
ously symmetric carbuncles is computed. The trigger mechanism remains the 
same: in a one-cell-high filament at the lower domain boundary, horizontal pre- 
shock velocity is reduced to 0. In similarity coordinates the resulting contact 
discontinuity is not preserved exactly; additional waves propagate upwards and 
the filament weakens noticeably (see Figure |3J|. However, this unelegant sideef- 
fect is not too large for our purposes. With the filament located at the lower 
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Figure 4. Filament emanating Figure 5. Filament and carbuncle tip (con- 
from the lowest edge in the left tact smeared due to poor grid alignment) 
boundary 




Figure 7. Detailed self-similar structure of carbuncle in Figure HO 




boundary, it has a lot of similarity with numerical boundary layers. Indeed, in 
the experience of the author these are the most frequent cause of carbuncle- like 
phenomena. 

In this new setting, carbuncle patterns can be computed in great detail (see 
Figure for overview and Figures [7| and [S] for more detail) . Better resolution 
could be achieved by using higher-order modifications, but here we choose the 
Godunov scheme because of its theoretical appeal (it manifestly satisfies all 
discrete entropy inequalities). 



4 Relation to entropy solution non-uniqueness 

In |EU06| a numerical example is presented that suggests that entropy solutions 
of the Cauchy problem for the compressible Euler equations are not always 
unique. More precisely, Figure [SJ constitutes an exact steady self-similar entropy 
solution; however, using it as initial data, numerical schemes produce another 
solution, Figure 1 101 which appears to be self-similar and (being computed by 
the Godunov scheme) an entropy solution, but clearly unsteady. 

The initial data in this example was motivated by research on the problem of 
supersonic flow onto a solid wedge. For small wedge opening angles (a in Figure 
EJ), there are two possible solutions, a strong and a (comparatively) weak shock 
(which was chosen for computing Figure HU|) . Replacing the solid wedge by gas 
with post-shock density and zero velocity yields Figure El 

In Figure [3 the tip of the carbuncle in Figure El is shown in detail. Clearly it 
resembles Figure El (the post-shock area in Figure El is supersonic which means 
the shock corresponds to the weak shock in the wedge problem). This anal- 
ogy suggests that carbuncles are another manifestation of the same problem: 
two different entropy solutions for identical data. In this case, the theoretical 
solution (the plane shock wave) is known to be physically correct as well as a 
vanishing viscosity and hence entropy solution; the carbuncle appears to be a 
second entropy solution which is probably unphysical without the presence of 
the filament. 

This paper is not the first to conjecture that carbuncles are phenomena of the 
continuum equations rather than purely numerical artefacts. Other references 
are |RGC.T0fi] where a linear stability analysis of plane shocks in the exact Euler 
equations seems to reveal a new unstable mode related to carbuncles f |CBGS02j 
have raised objections to this analysis), as well as |MRf)l| and |DMG03j both 
of which explicitly state suspicions that carbuncles and other anomalies are 
entropy solutions. 
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Figure 9. Probably unphysical steady self-similar entropy solution 




Figure 10. A second self-similar (but unsteady) entropy solution, for Figure El 
as initial data 
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It should be emphasized that, although the carbuncle in Figure El does not ap- 
pear to vanish as the grid is refined, it shrinks and grows and does not seem to 
converge to a single pattern either. Indeed, if the comparison to the nonunique- 
ness example is appropriate, then there may be a parametrized family of car- 
buncle patterns rather than a single one. There arc at least two parameters 
that are currently not under our control: the tip angle (3 and the location of 
the carbuncle tip. These are determined in some unknown and probably rather 
sensitive way by the numerical method, the filament and the inner structure 
of the carbuncle pattern. If one or both parameters are free (i.e. if there is a 
different carbuncle pattern for each choice of parameter in some interval), the 
carbuncle calculation should not be expected to converge in general. 

More generally, we should stress that even on a single fixed grid Figure [6] is 
rather unstable; small modifications to the numerical scheme can have drastic 
effects. For example, the Godunov scheme requires solving a nonlinear equation 
for every edge in every time step; solving this equation with less than machine 
accuracy can either trigger oscillations that abort the computation or produce 
additional smaller carbuncles along the shock that grow to spoil the overall 
pattern. 

It is possible to stabilize the calculation to some degree by adding numerical 
viscosity everywhere except near the lower boundary; this suppresses the ap- 
pearance of additional carbuncles, but smears the carbuncle pattern even more. 



5 Conclusions and open problems 

This article presents a new technique for triggering carbuncles in shock waves. 
Its main advantage is reliability; previous carbuncle computations were rather 
precarious and sensitive to changes in setting and numerical methods. There 
are two immediate benefits: first, it becomes clear that a rather large variety of 
numerical methods is affected, more than may have been suspected previously. 
Second, carbuncle structure can be computed in fine detail. This bears fruit 
immediately: as discussed in the previous section, Figure suggests a connec- 
tion between the carbuncle phenomenon and the hypothetical entropy solution 
nonuniqueness example (see above). Moreover, it can be verified that carbuncles 
are self-similar structures. Numerical experiments in which carbuncles appear 
to be steady rather than self-similar most likely exhibit equilibrum states in 
which further carbuncle growth would weaken numerical effects that stabilize 
the carbuncle tip. 

The increased understanding of the nature of carbuncles should facilitate the 
development of better criteria for predicting the appearance of carbuncles and 
of numerical methods that avoid carbuncles at a minimal loss of accuracy. Most 
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importantly, there is clear evidence that carbuncles are entropy solutions of the 
compressible Euler equations rather than non-entropy solutions or numerical 
artifacts. If carbuncle phenomena are indeed caused by a narrow family of con- 
tinuum functions, suppressing them might be easier than previously expected. 

Nevertheless, this work opens some new questions and leaves several unsolved 
problems. The most central ones arc closely related: first, it is not known how 
many parameters index the family of carbuncles. Drawing on the analogy to 
the nommiqueness example, there could be two parameters, tip angle (3 and 
tip location. It can be observed that these parameters change as the grid is 
refined which suggests that at least one of these parameters is free; there is no 
indication whether the other one is free or rigid. 

Moreover, it is desirable to have techniques for controlling the free parameters 
during computations. Using this control, the numerical scheme could be forced 
to converge to a single carbuncle. This would add substantial evidence for 
the hypothesis that carbuncles are entropy solutions: according to the Lax- 
Wendroff theorem (see |LW60I IKRW96I IGR961 tEllar) 1 . a convergent sequence 
of approximations generated by the Godunov scheme must have an entropy 
solution as limit. Moreover, with tight control of the free parameters, carbuncle 
structure could be computed faster, more precisely and more reliably. 

The internal structure in Figure[7]is quite smeared, so it will change significantly 
as the grid becomes finer (even in cases where the free parameters happen/are 
controlled to remain unchanged). Besides lack of control over free parame- 
ters, another obstacle towards computing the structure accurately is the strong 
smearing of contact discontinuities that are not grid-aligned. 

Even in the poor resolution of Figure[7|it is obvious that the carbuncle structure 
is rather complex. This might in part be due to the rather large pre-shock Mach 
number of 3.0 which was chosen for the pronounced carbuncle it causes. The 
carbuncle structure may be simpler for Mach numbers closer to 1, but for such 
values the author was not able to trigger a carbuncle in similarity coordinates. In 
standard coordinates, noticeable carbuncles start to appear at about Mach 1.4. 
It is not clear whether this number is somehow significant or merely a limitation 
of the filament trigger mechanism. The carbuncles grow very slowly near Mach 
1.4, so they may simply escape observation for smaller Mach numbers. 

In Figures |5] and 1101 there appear to be two different entropy solutions. This 
by itself does not rule out that there are only a few isolated entropy solutions; 
in this case the entropy condition would still be rather useful. However, if the 
tip shock opening angle is a free parameter that can be shifted up to 180° (this 
requires a strong shock at the carbuncle tip and causes the angle a between 
the contacts to decrease to 0°), there may be a continuous (in L 1 ) family of 
carbuncle patterns, containing the physically correct plane shock wave as a limit 
case — such a continuum of entropy solutions for the same data would seriously 
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reduce the value of the entropy condition. Numerical calculations can be used 
to verify these ideas once we find techniques to control the tip angle j3 (this is 
another motivation for trying to control the free parameters). However, due to 
this theoretical relevance of carbuncles, it would be desirable to prove existence 
of such patterns mathematically rather than numerically. Due to the complex 
structure of the carbuncle patterns and the intricacies of the Euler equations, 
this appears to be extremely difficult. 



References 

[Ais04] H. Aiso, Instability analysis in conservative difference approxima- 
tions for compressible euler equations, Presentation at 10th Int'l 
Conference on Hyperbolic Problems Theory, Numerics, Applications 
(HYP2004), Osaka, Japan, 2004. 

[CBGS02] J.-F. Coulombel, S. Benzoni-Gavage, and D. Serre, Note on a pa- 
per by Robinet, Gressier, Casalis & Moschetta, J. Fluid Mech. 469 
(2002), 401-405. 

[DMG03] M. Dumbser, J.-M. Moschetta, and J. Gressier, A matrix stability 
analysis of the carbuncle phenomenon, submitted to Elsevier Science, 
2003. 

[E1100] V. Elling, Numerical simulation of gas flow in moving domains, 
Diploma Thesis, RWTH Aachen (Germany), 2000, http://www- 
sccm.stanford.edu/~clling/diplom-abstract.html. 

[E1105] , A Lax- Wendroff type theorem for unstructured grids, 

Ph.D. Dissertation, Stanford University, 2005, http://www- 
sccm.stanford.edu/~elling/phd-abstract.html. 

[E1106] , A possible counterexample to well-posedness of entropy so- 
lution and to Godunov scheme convergence, Math. Comp. 75 (2006), 
1721-1733, see also arxiv:math.NA/0509331. 

[Ellar] , A Lax- Wendroff type theorem for unstructured 

quasi-uniform grids, Math. Comp. (to appear), see also 
arxiv:math.NA/0509333. 

[God59] S. K. Godunov, A finite difference method for the numerical compu- 
tation of discontinuous solutions of the equations of fluid dynamics, 
Mat. Sb. 47 (1959), 271-290. 

[GR96] E. Godlewski and P.-A. Raviart, Numerical approximation of hyper- 
bolic systems of conservation laws, Springer, 1996. 



11 



[KRW96] D. Kroner, M. Rokyta, and M. Wierse, A Lax- Wendroff type theorem 
for upwind finite volume schemes in 2-D, East- West J. Numer. Math. 
4 (1996), 279-292. 

[LW60] P. Lax and B. Wendroff, Systems of conservation laws, Comm. Pure 
Appl. Math. 13 (1960), 217-237. 

[MR01] K.W. Morton and P. Roe, Vorticity-preserving Lax-Wendroff-type 
schemes for the system wave equation, SIAM J. Sci. Comput. 23 
(2001), no. 1, 170-192. 

[OS82] S. Osher and F. Solomon, Upwind difference schemes for hyperbolic 
systems of conservation laws, Math. Comp. 38 (1982), 339-373. 

[PI88] KM. Peery and S.T. Imlay, Blunt-body flow simulations, AIAA paper 
88-2904, 1988. 

[Qui94] J. Quirk, A contribution to the great Riemann solver debate, Intl. J. 
Numer. Mcth. Fluids 18 (1994), 555-574. 

[RGCJ00] J.-Ch. Robinet, J. Gressier, G. Casalis, and J.-M. Moschetta, Shock 
wave instability and the carbuncle phenomenon: same intrinsic ori- 
gin?, J. Fluid Mech. 417 (2000), 237-263. 

[Roe81] P.L. Roc, Approximate Riemann solvers, parameter vectors, and dif- 
ference schemes, J. Comput. Phys. 43 (1981), 357-372. 

[S089] C. W. Shu and S. Osher, Efficient implementation of essentially non- 
oscillatory shock- capturing schemes, II, J. Comp. Phys. 83 (1989), 
32-78. 



12 



